Tomato for Carlos
library(rstan)
library(rethinking)
library(brms)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
if(file.exists("~/Google Drive/Rdata/carlos.Rdata")) load("~/Google Drive/Rdata/carlos.Rdata")tomato <- read.csv("TomatoR2CSHL.csv")
summary(tomato) shelf flat col row acs trt days date
U:161 Min. : 1.00 G :133 Min. :1.00 LA1954 : 40 H:495 Min. :28.00 5/5/08:716
V:174 1st Qu.: 9.00 H :127 1st Qu.:2.00 LA2695 : 39 L:513 1st Qu.:28.00 5/6/08:292
W:178 Median :17.00 F :125 Median :3.00 LA1361 : 37 Median :28.00
X:174 Mean :17.89 C :117 Mean :2.56 LA2167 : 37 Mean :28.29
Y:125 3rd Qu.:28.00 D :117 3rd Qu.:4.00 LA2773 : 37 3rd Qu.:29.00
Z:196 Max. :36.00 E :107 Max. :4.00 LA1474 : 36 Max. :29.00
(Other):282 (Other):782
hyp int1 int2 int3 int4 intleng
Min. : 6.17 Min. : 0.00 Min. : 0.000 Min. : 0.010 Min. : 0.030 Min. : 0.000
1st Qu.:26.81 1st Qu.: 1.74 1st Qu.: 1.060 1st Qu.: 2.975 1st Qu.: 2.163 1st Qu.: 9.637
Median :32.02 Median : 3.59 Median : 3.120 Median : 5.625 Median : 3.995 Median :17.255
Mean :33.36 Mean : 4.71 Mean : 4.287 Mean : 6.794 Mean : 5.102 Mean :20.340
3rd Qu.:38.56 3rd Qu.: 6.46 3rd Qu.: 6.320 3rd Qu.: 9.367 3rd Qu.: 7.018 3rd Qu.:28.145
Max. :74.60 Max. :39.01 Max. :28.980 Max. :27.760 Max. :23.280 Max. :92.420
NA's :1 NA's :1 NA's :4 NA's :102
totleng petleng leafleng leafwid leafnum ndvi
Min. : 13.59 Min. : 1.53 Min. : 9.74 Min. : 8.29 Min. :3.000 Min. :100
1st Qu.: 39.25 1st Qu.:11.20 1st Qu.:27.43 1st Qu.:29.48 1st Qu.:5.000 1st Qu.:108
Median : 50.98 Median :15.13 Median :34.59 Median :39.62 Median :5.000 Median :115
Mean : 53.70 Mean :15.92 Mean :35.54 Mean :39.29 Mean :5.063 Mean :118
3rd Qu.: 64.76 3rd Qu.:20.48 3rd Qu.:42.98 3rd Qu.:47.75 3rd Qu.:6.000 3rd Qu.:128
Max. :129.43 Max. :44.44 Max. :95.19 Max. :90.27 Max. :8.000 Max. :137
NA's :2 NA's :1 NA's :1 NA's :1
lat lon alt species who
Min. :-25.400 Min. :-78.52 Min. : 0 S. chilense :207 Dan :402
1st Qu.:-16.607 1st Qu.:-75.92 1st Qu.:1020 S. chmielewskii:226 Pepe:606
Median :-14.152 Median :-73.63 Median :2240 S. habrochaites:226
Mean :-14.490 Mean :-73.71 Mean :2035 S. pennellii :132
3rd Qu.:-12.450 3rd Qu.:-71.70 3rd Qu.:3110 S. peruvianum :217
Max. : -5.767 Max. :-68.07 Max. :3540
hyp1 <- brm(hyp ~ acs*trt + (trt|species) + (1|shelf),
prior = c(
set_prior("normal(20,10)",class="Intercept"),
set_prior("normal(0,10)",class="b"),
set_prior("cauchy(0,1)",class="sigma")
),
iter=5000,
warmup=1500,
data=tomato)save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp2 <- brm(hyp ~ acs*trt + (1|species) + (1|shelf),
prior = c(
set_prior("normal(20,10)",class="Intercept"),
set_prior("normal(0,10)",class="b"),
set_prior("cauchy(0,1)",class="sigma")
),
iter=5000,
warmup=1500,
data=tomato)Compiling the C++ model
Start samplingstarting worker pid=99011 on localhost:11828 at 12:22:36.419
starting worker pid=99019 on localhost:11828 at 12:22:36.592
starting worker pid=99027 on localhost:11828 at 12:22:36.760
starting worker pid=99035 on localhost:11828 at 12:22:36.935
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 1501 / 5000 [ 30%] (Sampling)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 1501 / 5000 [ 30%] (Sampling)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1501 / 5000 [ 30%] (Sampling)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2, Iteration: 1501 / 5000 [ 30%] (Sampling)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Sampling)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 180.896 seconds (Warm-up)
226.983 seconds (Sampling)
407.879 seconds (Total)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 180.073 seconds (Warm-up)
230.852 seconds (Sampling)
410.924 seconds (Total)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 208.94 seconds (Warm-up)
243.667 seconds (Sampling)
452.607 seconds (Total)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 171.104 seconds (Warm-up)
349.297 seconds (Sampling)
520.401 seconds (Total)There were 541 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp3 <- update(hyp1, formula. = hyp ~ acs*trt + (1|species))The desired formula changes require recompling the model
Compiling the C++ model
Start samplingstarting worker pid=99334 on localhost:11828 at 12:32:57.610
starting worker pid=99342 on localhost:11828 at 12:32:57.791
starting worker pid=99350 on localhost:11828 at 12:32:57.969
starting worker pid=99358 on localhost:11828 at 12:32:58.153
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 195.425 seconds (Warm-up)
63.3163 seconds (Sampling)
258.741 seconds (Total)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 199.978 seconds (Warm-up)
107.182 seconds (Sampling)
307.16 seconds (Total)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 205.032 seconds (Warm-up)
114.488 seconds (Sampling)
319.52 seconds (Total)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 207.932 seconds (Warm-up)
119.223 seconds (Sampling)
327.155 seconds (Total)There were 1105 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp4 <- update(hyp1, formula. = hyp ~ acs*trt + (1|shelf))The desired formula changes require recompling the model
Compiling the C++ model
Start samplingstarting worker pid=99586 on localhost:11828 at 12:39:50.457
starting worker pid=99594 on localhost:11828 at 12:39:50.613
starting worker pid=99602 on localhost:11828 at 12:39:50.770
starting worker pid=99610 on localhost:11828 at 12:39:50.922
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 151.259 seconds (Warm-up)
78.1472 seconds (Sampling)
229.407 seconds (Total)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 138.363 seconds (Warm-up)
105.666 seconds (Sampling)
244.029 seconds (Total)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 153.591 seconds (Warm-up)
94.0214 seconds (Sampling)
247.613 seconds (Total)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 144.983 seconds (Warm-up)
136.555 seconds (Sampling)
281.538 seconds (Total)There were 7 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp5 <- update(hyp1, formula. = hyp ~ acs*trt)The desired formula changes require recompling the model
Compiling the C++ model
Start samplingstarting worker pid=99715 on localhost:11828 at 12:45:40.313
starting worker pid=99723 on localhost:11828 at 12:45:40.486
starting worker pid=99731 on localhost:11828 at 12:45:40.654
starting worker pid=99739 on localhost:11828 at 12:45:40.814
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 85.807 seconds (Warm-up)
39.2256 seconds (Sampling)
125.033 seconds (Total)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 86.8877 seconds (Warm-up)
36.1756 seconds (Sampling)
123.063 seconds (Total)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 91.221 seconds (Warm-up)
36.5213 seconds (Sampling)
127.742 seconds (Total)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 90.6127 seconds (Warm-up)
38.1106 seconds (Sampling)
128.723 seconds (Total)Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")loo(hyp1,hyp2,hyp3,hyp4,hyp5)Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details. LOOIC SE
hyp1 7042.25 58.38
hyp2 7042.55 58.24
hyp3 7088.25 58.10
hyp4 7041.43 58.26
hyp5 7084.61 58.15
hyp1 - hyp2 -0.30 1.00
hyp1 - hyp3 -46.00 12.92
hyp1 - hyp4 0.82 1.07
hyp1 - hyp5 -42.36 12.85
hyp2 - hyp3 -45.70 12.78
hyp2 - hyp4 1.11 0.48
hyp2 - hyp5 -42.06 12.71
hyp3 - hyp4 46.81 12.76
hyp3 - hyp5 3.64 0.84
hyp4 - hyp5 -43.18 12.68So we should go with hyp4 (no species effect)
summary(hyp4) Family: gaussian (identity)
Formula: hyp ~ acs + trt + (1 | shelf) + acs:trt
Data: tomato (Number of observations: 1008)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
WAIC: Not computed
Group-Level Effects:
~shelf (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 2.85 1.34 1.23 6.36 3387 1
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 28.83 2.09 24.57 32.97 3863 1
acsLA1305 -0.05 2.45 -4.96 4.65 4684 1
acsLA1306 -1.56 2.22 -5.91 2.82 5142 1
acsLA1316 9.38 2.24 5.01 13.74 4346 1
acsLA1317 5.73 2.16 1.51 9.92 4888 1
acsLA1318 -1.88 2.19 -6.23 2.47 4024 1
acsLA1336 11.36 2.29 6.81 15.83 4572 1
acsLA1361 1.46 2.07 -2.76 5.50 4150 1
acsLA1363 -0.67 2.16 -4.92 3.57 4379 1
acsLA1367 -8.24 2.81 -13.80 -2.75 10000 1
acsLA1474 12.64 2.13 8.49 16.80 4482 1
acsLA1691 -1.79 2.31 -6.34 2.76 4669 1
acsLA1912 -5.45 2.57 -10.55 -0.45 5795 1
acsLA1918 6.63 2.30 2.15 11.22 5424 1
acsLA1928 0.71 3.97 -6.95 8.65 10000 1
acsLA1954 7.32 2.05 3.34 11.33 3930 1
acsLA1958 7.44 2.16 3.24 11.69 4442 1
acsLA1969 -4.93 2.26 -9.37 -0.49 4521 1
acsLA1973 6.00 2.32 1.48 10.42 5604 1
acsLA2167 -1.29 2.08 -5.38 2.78 3931 1
acsLA2204 -1.73 2.05 -5.76 2.24 3863 1
acsLA2560 -5.71 3.88 -13.30 1.91 10000 1
acsLA2580 -6.53 2.67 -11.78 -1.28 5300 1
acsLA2663 -2.22 2.14 -6.40 1.96 4201 1
acsLA2695 0.31 2.08 -3.76 4.37 4053 1
acsLA2748 -1.61 2.56 -6.60 3.37 5837 1
acsLA2773 4.64 2.09 0.59 8.79 4160 1
acsLA2884 1.25 2.34 -3.31 5.81 4433 1
acsLA2891 0.65 2.38 -4.08 5.24 5447 1
acsLA2931 8.31 2.30 3.81 12.82 4775 1
acsLA361 2.72 2.36 -1.85 7.34 4959 1
acsLA3788 -4.77 2.27 -9.19 -0.32 4992 1
acsLA3791 -1.46 2.45 -6.24 3.32 5022 1
acsLA3795 0.67 2.24 -3.73 4.98 4991 1
acsLA446 1.57 2.36 -3.08 6.17 4771 1
acsLA716 -9.74 3.46 -16.41 -2.87 10000 1
trtL 3.66 2.76 -2.09 9.09 3947 1
acsLA1305:trtL -5.66 3.16 -11.85 0.45 6853 1
acsLA1306:trtL -1.23 2.86 -6.94 4.36 6470 1
acsLA1316:trtL 0.88 2.90 -4.76 6.48 6674 1
acsLA1317:trtL -1.38 2.82 -6.92 4.17 10000 1
acsLA1318:trtL -0.44 2.95 -6.16 5.45 6977 1
acsLA1336:trtL 3.79 2.93 -1.95 9.53 6621 1
acsLA1361:trtL 2.74 2.78 -2.77 8.22 10000 1
acsLA1363:trtL -2.44 2.86 -7.97 3.25 6157 1
acsLA1367:trtL 9.39 3.78 2.00 16.80 10000 1
acsLA1474:trtL 3.83 2.80 -1.69 9.15 10000 1
acsLA1691:trtL -0.61 3.14 -6.77 5.48 10000 1
acsLA1912:trtL 0.30 3.26 -6.07 6.52 6838 1
acsLA1918:trtL -1.63 3.26 -8.04 4.82 10000 1
acsLA1928:trtL -9.85 5.75 -21.04 1.63 10000 1
acsLA1954:trtL 1.68 2.68 -3.58 6.96 6278 1
acsLA1958:trtL 2.17 2.89 -3.56 7.81 10000 1
acsLA1969:trtL -0.88 3.08 -6.85 5.11 5883 1
acsLA1973:trtL 8.11 3.02 2.29 14.08 10000 1
acsLA2167:trtL -0.54 2.78 -6.07 4.96 5263 1
acsLA2204:trtL 1.20 2.83 -4.48 6.68 10000 1
acsLA2560:trtL 6.76 4.80 -2.62 16.17 10000 1
acsLA2580:trtL 4.17 3.38 -2.30 10.85 10000 1
acsLA2663:trtL -1.37 2.81 -6.96 4.08 5960 1
acsLA2695:trtL -0.78 2.74 -6.11 4.60 6219 1
acsLA2748:trtL -0.40 3.30 -6.95 6.13 10000 1
acsLA2773:trtL 0.09 2.80 -5.54 5.44 6499 1
acsLA2884:trtL 1.70 3.07 -4.40 7.64 10000 1
acsLA2891:trtL 7.06 3.20 0.91 13.22 10000 1
acsLA2931:trtL 4.77 2.94 -0.98 10.42 10000 1
acsLA361:trtL 0.80 3.05 -5.03 6.76 10000 1
acsLA3788:trtL 10.09 3.21 3.73 16.45 10000 1
acsLA3791:trtL 3.37 3.09 -2.80 9.39 10000 1
acsLA3795:trtL 3.31 2.95 -2.51 9.12 10000 1
acsLA446:trtL -0.47 3.27 -6.88 5.88 10000 1
acsLA716:trtL 4.76 4.35 -3.81 13.07 10000 1
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 7.67 0.18 7.34 8.03 10000 1
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(hyp4, ask=FALSE)
ranef(hyp4)$shelf
Intercept
U -1.4409542
V 1.7748644
W 0.6353095
X -2.0373836
Y -0.1921667
Z 2.4487958pred.df <- unique(tomato[,c("species","acs","trt")])
head(pred.df)hyp.prediction <- cbind(pred.df,fitted(hyp4,newdata = pred.df, re_formula = NA))
colnames(hyp.prediction)[6:7] <- c("low95","high95")
head(hyp.prediction)write.csv(hyp.prediction,file="hyp_accession.csv")pl <- ggplot(hyp.prediction,aes(x=acs,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl <- pl + facet_wrap( ~ species, ncol=2, scales="free_x")
pl + ggtitle("hypocotyls per accession")
ggsave("hypocotyls_accession.pdf",height=11,width=8.5)
by(tomato$intleng,tomato$trt,mean)
int1 <- brm(intleng ~ acs*trt + (trt|species) + (1|shelf),
prior = c(
set_prior("normal(10,10)",class="Intercept"),
set_prior("normal(0,20)",class="b"),
set_prior("cauchy(0,1)",class="sigma")
),
iter=5000,
warmup=1500,
data=tomato)save.image(file="~/Google Drive/Rdata/carlos.Rdata")int2 <- update(int1, formula. = intleng ~ acs*trt + (1|species) + (1|shelf))The desired formula changes require recompling the model
Compiling the C++ model
Start samplingstarting worker pid=177 on localhost:11828 at 13:03:31.512
starting worker pid=185 on localhost:11828 at 13:03:31.684
starting worker pid=197 on localhost:11828 at 13:03:31.854
starting worker pid=205 on localhost:11828 at 13:03:32.026
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 228.748 seconds (Warm-up)
165.635 seconds (Sampling)
394.382 seconds (Total)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 238.932 seconds (Warm-up)
163.406 seconds (Sampling)
402.338 seconds (Total)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 258.535 seconds (Warm-up)
149.635 seconds (Sampling)
408.17 seconds (Total)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 229.139 seconds (Warm-up)
265.832 seconds (Sampling)
494.97 seconds (Total)There were 533 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")int3 <- update(int2, formula. = intleng ~ acs*trt + (1|species))The desired formula changes require recompling the model
Compiling the C++ model
Start samplingstarting worker pid=502 on localhost:11828 at 13:16:43.400
starting worker pid=515 on localhost:11828 at 13:16:43.568
starting worker pid=524 on localhost:11828 at 13:16:43.724
starting worker pid=532 on localhost:11828 at 13:16:43.880
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 222.004 seconds (Warm-up)
31.0146 seconds (Sampling)
253.019 seconds (Total)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 206.206 seconds (Warm-up)
76.8814 seconds (Sampling)
283.087 seconds (Total)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 212.625 seconds (Warm-up)
117.46 seconds (Sampling)
330.086 seconds (Total)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 240.151 seconds (Warm-up)
95.9457 seconds (Sampling)
336.097 seconds (Total)There were 2407 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")int4 <- update(int3, formula. = intleng ~ acs*trt + (1|shelf))The desired formula changes require recompling the model
Compiling the C++ model
Start samplingstarting worker pid=838 on localhost:11828 at 13:27:07.734
starting worker pid=848 on localhost:11828 at 13:27:07.915
starting worker pid=863 on localhost:11828 at 13:27:08.097
starting worker pid=875 on localhost:11828 at 13:27:08.272
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 126.289 seconds (Warm-up)
70.5382 seconds (Sampling)
196.827 seconds (Total)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 133.997 seconds (Warm-up)
71.3704 seconds (Sampling)
205.367 seconds (Total)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 129.769 seconds (Warm-up)
79.0085 seconds (Sampling)
208.778 seconds (Total)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 138.738 seconds (Warm-up)
78.9029 seconds (Sampling)
217.641 seconds (Total)There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")int5 <- update(int1, formula. = intleng ~ acs*trt)The desired formula changes require recompling the model
Compiling the C++ model
Start samplingstarting worker pid=1172 on localhost:11828 at 13:32:53.126
starting worker pid=1180 on localhost:11828 at 13:32:53.309
starting worker pid=1188 on localhost:11828 at 13:32:53.487
starting worker pid=1196 on localhost:11828 at 13:32:53.668
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4, Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 90.9246 seconds (Warm-up)
44.8574 seconds (Sampling)
135.782 seconds (Total)
Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 97.5076 seconds (Warm-up)
49.4388 seconds (Sampling)
146.946 seconds (Total)
Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 96.0257 seconds (Warm-up)
62.106 seconds (Sampling)
158.132 seconds (Total)
Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
Elapsed Time: 86.7937 seconds (Warm-up)
75.9581 seconds (Sampling)
162.752 seconds (Total)Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2
2: package ‘ggplot2’ was built under R version 3.3.2
3: package ‘StanHeaders’ was built under R version 3.3.2 save.image(file="~/Google Drive/Rdata/carlos.Rdata")loo(int1,int2,int3,int4,int5)Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details. LOOIC SE
int1 7396.58 64.57
int2 7399.29 64.61
int3 7413.85 62.36
int4 7397.58 64.51
int5 7412.23 62.83
int1 - int2 -2.71 0.92
int1 - int3 -17.28 10.76
int1 - int4 -1.00 0.82
int1 - int5 -15.66 10.50
int2 - int3 -14.56 10.70
int2 - int4 1.71 0.65
int2 - int5 -12.94 10.44
int3 - int4 16.28 10.59
int3 - int5 1.62 2.20
int4 - int5 -14.66 10.33So we should go with int4 (no species effect). Case could also be made for int5 (no random effects)
summary(int4) Family: gaussian (identity)
Formula: intleng ~ acs + trt + (1 | shelf) + acs:trt
Data: tomato (Number of observations: 1008)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
WAIC: Not computed
Group-Level Effects:
~shelf (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 2.35 1.33 0.82 5.72 2915 1
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 8.80 2.48 3.79 13.60 1058 1
acsLA1305 3.38 3.29 -2.98 9.89 1902 1
acsLA1306 1.99 3.02 -3.92 7.70 1618 1
acsLA1316 0.79 3.04 -5.17 6.75 1603 1
acsLA1317 11.85 3.00 5.97 17.75 1708 1
acsLA1318 7.79 3.00 1.83 13.57 1649 1
acsLA1336 7.23 3.10 1.08 13.32 1678 1
acsLA1361 5.39 2.83 -0.17 10.96 1542 1
acsLA1363 11.02 2.93 5.14 16.80 1593 1
acsLA1367 0.15 3.65 -6.99 7.42 2263 1
acsLA1474 1.50 2.87 -4.21 7.08 1518 1
acsLA1691 1.20 3.10 -4.85 7.20 1665 1
acsLA1912 -2.30 3.47 -9.10 4.49 2375 1
acsLA1918 3.29 3.14 -2.92 9.55 1739 1
acsLA1928 2.26 5.24 -7.92 12.45 4792 1
acsLA1954 6.35 2.82 0.83 11.87 1457 1
acsLA1958 -1.35 2.96 -7.15 4.47 1670 1
acsLA1969 -6.59 3.06 -12.63 -0.54 1628 1
acsLA1973 11.55 3.10 5.50 17.64 1627 1
acsLA2167 3.58 2.87 -2.08 9.14 1573 1
acsLA2204 8.46 2.83 2.92 13.96 1368 1
acsLA2560 0.13 5.20 -10.02 10.33 4298 1
acsLA2580 -4.68 3.56 -11.68 2.28 2107 1
acsLA2663 4.95 2.96 -0.88 10.74 1573 1
acsLA2695 3.41 2.87 -2.33 8.95 1427 1
acsLA2748 -2.50 3.48 -9.50 4.35 2026 1
acsLA2773 -3.84 2.87 -9.61 1.74 1532 1
acsLA2884 -2.37 3.08 -8.39 3.65 1658 1
acsLA2891 -2.03 3.18 -8.30 4.17 1804 1
acsLA2931 -3.28 3.08 -9.27 2.79 1627 1
acsLA361 12.37 3.20 5.92 18.58 1791 1
acsLA3788 -3.55 3.10 -9.65 2.65 1627 1
acsLA3791 -6.13 3.27 -12.57 0.29 1837 1
acsLA3795 -0.67 3.04 -6.65 5.26 1604 1
acsLA446 2.15 3.11 -4.02 8.18 1657 1
acsLA716 -2.32 4.71 -11.61 6.92 3806 1
trtL 12.68 3.25 6.06 19.04 1177 1
acsLA1305:trtL 1.89 4.22 -6.43 10.14 2225 1
acsLA1306:trtL 5.07 3.90 -2.55 12.83 1968 1
acsLA1316:trtL -1.84 3.96 -9.49 6.02 1867 1
acsLA1317:trtL 7.93 3.95 0.40 15.77 1736 1
acsLA1318:trtL 6.87 4.05 -1.05 14.91 2048 1
acsLA1336:trtL 7.23 4.07 -0.80 15.20 1797 1
acsLA1361:trtL 15.69 3.85 8.21 23.35 1884 1
acsLA1363:trtL 5.65 3.90 -1.89 13.30 1866 1
acsLA1367:trtL 8.91 4.97 -0.95 18.60 3153 1
acsLA1474:trtL 0.51 3.82 -6.87 8.09 1778 1
acsLA1691:trtL -5.80 4.22 -14.15 2.39 2193 1
acsLA1912:trtL 2.23 4.44 -6.38 10.99 2420 1
acsLA1918:trtL -5.97 4.44 -14.55 2.65 2279 1
acsLA1928:trtL -1.75 7.81 -17.15 13.58 10000 1
acsLA1954:trtL 1.84 3.72 -5.36 9.26 1746 1
acsLA1958:trtL 6.58 3.96 -1.24 14.34 2045 1
acsLA1969:trtL 6.55 4.15 -1.56 14.84 2226 1
acsLA1973:trtL 15.81 4.05 7.95 23.79 2075 1
acsLA2167:trtL 0.88 3.84 -6.64 8.41 1815 1
acsLA2204:trtL 5.15 3.90 -2.41 12.78 1836 1
acsLA2560:trtL 10.77 6.49 -1.93 23.36 4647 1
acsLA2580:trtL 8.10 4.60 -0.83 17.04 2313 1
acsLA2663:trtL -0.65 3.90 -8.23 7.11 1835 1
acsLA2695:trtL -4.27 3.76 -11.62 3.19 1683 1
acsLA2748:trtL 17.95 4.54 8.96 26.81 2430 1
acsLA2773:trtL 13.28 3.83 5.76 20.82 1765 1
acsLA2884:trtL 0.50 4.09 -7.45 8.55 1936 1
acsLA2891:trtL 8.71 4.23 0.34 16.92 2154 1
acsLA2931:trtL 9.45 3.96 1.67 17.27 1823 1
acsLA361:trtL 8.33 4.20 0.15 16.69 2112 1
acsLA3788:trtL 4.89 4.29 -3.43 13.33 2042 1
acsLA3791:trtL 1.49 4.15 -6.51 9.71 2233 1
acsLA3795:trtL 4.81 4.02 -3.23 12.67 1841 1
acsLA446:trtL 6.34 4.32 -2.07 14.82 2257 1
acsLA716:trtL -0.28 6.05 -12.14 11.42 3960 1
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 9.12 0.21 8.72 9.54 10000 1
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(int4, ask=FALSE)
ranef(int4)$shelf
Intercept
U -1.60562115
V 2.34342216
W -0.27581150
X -0.01534482
Y -0.56118509
Z 0.73463441intleng.prediction <- cbind(pred.df,fitted(int4,newdata = pred.df, re_formula = NA))
colnames(intleng.prediction)[6:7] <- c("low95","high95")
head(intleng.prediction)write.csv(intleng.prediction,file="intleng_accession.csv")pl <- ggplot(intleng.prediction,aes(x=acs,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl <- pl + facet_wrap( ~ species, ncol=2, scales="free_x")
pl + ggtitle("internode")
ggsave("internodes_accession.pdf",height=11,width=8.5)
species level analyses
hypocotyl
hyp6 <- brm(hyp ~ species*trt + (trt|acs) + (1|shelf),
prior = c(
set_prior("normal(20,10)",class="Intercept"),
set_prior("normal(0,10)",class="b"),
set_prior("cauchy(0,1)",class="sigma")
),
iter=5000,
warmup=1500,
data=tomato)save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp7 <- update(hyp6, formula. = . ~ species*trt + (1|acs) + (1|shelf) )save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp8 <- update(hyp7, formula. = . ~ species*trt + (1|shelf) )save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp9 <- update(hyp7, formula. = . ~ species*trt + (1|acs) )save.image(file="~/Google Drive/Rdata/carlos.Rdata")hyp10 <- update(hyp7, formula. = . ~ species*trt)save.image(file="~/Google Drive/Rdata/carlos.Rdata")summary(hyp6) Family: gaussian (identity)
Formula: hyp ~ species * trt + (trt | acs) + (1 | shelf)
Data: tomato (Number of observations: 1008)
Samples: 4 chains, each with iter = 5000; warmup = 1500; thin = 1;
total post-warmup samples = 14000
WAIC: Not computed
Group-Level Effects:
~acs (Number of levels: 36)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 4.15 0.69 2.99 5.68 5208 1
sd(trtL) 1.81 0.71 0.45 3.28 4807 1
cor(Intercept,trtL) 0.67 0.27 0.01 0.99 6552 1
~shelf (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 2.88 1.37 1.24 6.49 4139 1
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 30.90 2.46 25.93 35.65 3643 1
speciesS.chmielewskii -0.51 2.38 -5.14 4.19 4060 1
speciesS.habrochaites -1.60 2.32 -6.17 3.01 3681 1
speciesS.pennellii -8.41 2.49 -13.26 -3.47 4310 1
speciesS.peruvianum 3.52 2.40 -1.19 8.24 3671 1
trtL 5.50 2.68 -0.07 10.59 4802 1
speciesS.chmielewskii:trtL -3.04 1.72 -6.44 0.33 8131 1
speciesS.habrochaites:trtL -2.30 1.71 -5.72 1.09 8232 1
speciesS.pennellii:trtL 4.31 1.93 0.46 8.12 8650 1
speciesS.peruvianum:trtL 0.19 1.73 -3.19 3.60 7729 1
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 7.67 0.18 7.33 8.02 14000 1
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).summary(hyp7) Family: gaussian (identity)
Formula: hyp ~ species + trt + (1 | acs) + (1 | shelf) + species:trt
Data: tomato (Number of observations: 1008)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
WAIC: Not computed
Group-Level Effects:
~acs (Number of levels: 36)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 4.85 0.69 3.69 6.4 2597 1
~shelf (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 3.02 1.7 1.22 7.53 627 1.01
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 30.83 2.54 25.74 35.77 3164 1
speciesS.chmielewskii -0.48 2.63 -5.66 4.78 3295 1
speciesS.habrochaites -1.44 2.50 -6.35 3.33 3365 1
speciesS.pennellii -8.34 2.66 -13.58 -3.04 3809 1
speciesS.peruvianum 3.49 2.57 -1.56 8.53 3504 1
trtL 5.56 2.74 -0.34 10.83 1666 1
speciesS.chmielewskii:trtL -3.19 1.45 -6.02 -0.33 6929 1
speciesS.habrochaites:trtL -2.38 1.45 -5.23 0.45 4668 1
speciesS.pennellii:trtL 4.17 1.64 0.99 7.44 6833 1
speciesS.peruvianum:trtL 0.26 1.44 -2.63 3.04 7030 1
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 7.71 0.18 7.37 8.06 10000 1
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).summary(hyp8) Family: gaussian (identity)
Formula: hyp ~ species + trt + (1 | shelf) + species:trt
Data: tomato (Number of observations: 1008)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
WAIC: Not computed
Group-Level Effects:
~shelf (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 2.53 1.23 1.01 5.78 2411 1
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 31.47 1.70 27.97 34.90 3964 1
speciesS.chmielewskii -1.06 1.18 -3.38 1.26 6018 1
speciesS.habrochaites -2.06 1.14 -4.29 0.13 5576 1
speciesS.pennellii -8.47 1.38 -11.20 -5.81 6226 1
speciesS.peruvianum 3.66 1.18 1.32 5.98 5608 1
trtL 6.01 2.40 1.11 10.73 3406 1
speciesS.chmielewskii:trtL -3.23 1.63 -6.44 -0.04 5875 1
speciesS.habrochaites:trtL -2.67 1.62 -5.85 0.45 5471 1
speciesS.pennellii:trtL 3.49 1.86 -0.10 7.16 6458 1
speciesS.peruvianum:trtL 0.03 1.63 -3.14 3.22 5323 1
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 8.88 0.2 8.5 9.27 10000 1
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).summary(hyp9) Family: gaussian (identity)
Formula: hyp ~ species + trt + (1 | acs) + species:trt
Data: tomato (Number of observations: 1008)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
WAIC: Not computed
Group-Level Effects:
~acs (Number of levels: 36)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 4.74 0.69 3.6 6.29 2491 1
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 30.98 1.81 27.39 34.58 1700 1
speciesS.chmielewskii -0.44 2.57 -5.54 4.59 2024 1
speciesS.habrochaites -1.50 2.50 -6.40 3.47 1920 1
speciesS.pennellii -7.64 2.68 -13.02 -2.38 2039 1
speciesS.peruvianum 3.51 2.59 -1.68 8.53 1903 1
trtL 6.07 1.05 4.02 8.13 2498 1
speciesS.chmielewskii:trtL -3.33 1.48 -6.22 -0.47 3381 1
speciesS.habrochaites:trtL -2.62 1.46 -5.49 0.23 3241 1
speciesS.pennellii:trtL 3.58 1.71 0.23 6.91 4482 1
speciesS.peruvianum:trtL 0.20 1.49 -2.72 3.07 3708 1
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 7.9 0.18 7.56 8.26 9577 1
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).summary(hyp10) Family: gaussian (identity)
Formula: hyp ~ species + trt + species:trt
Data: tomato (Number of observations: 1008)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
WAIC: Not computed
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 31.53 0.82 29.92 33.15 4801 1
speciesS.chmielewskii -0.98 1.15 -3.20 1.27 5388 1
speciesS.habrochaites -2.07 1.12 -4.28 0.12 5441 1
speciesS.pennellii -7.82 1.34 -10.42 -5.20 5804 1
speciesS.peruvianum 3.72 1.17 1.37 6.03 5668 1
trtL 6.28 1.14 4.05 8.53 4160 1
speciesS.chmielewskii:trtL -3.32 1.60 -6.48 -0.24 4773 1
speciesS.habrochaites:trtL -2.86 1.59 -5.93 0.23 4837 1
speciesS.pennellii:trtL 3.08 1.83 -0.53 6.66 5175 1
speciesS.peruvianum:trtL 0.00 1.62 -3.16 3.19 5134 1
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 9 0.19 8.63 9.39 8766 1
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).loo(hyp6,hyp7,hyp8,hyp9,hyp10) LOOIC SE
hyp6 7019.05 58.82
hyp7 7022.09 58.24
hyp8 7277.38 58.08
hyp9 7068.53 58.00
hyp10 7300.36 57.60
hyp6 - hyp7 -3.04 6.03
hyp6 - hyp8 -258.33 31.40
hyp6 - hyp9 -49.48 14.03
hyp6 - hyp10 -281.31 32.13
hyp7 - hyp8 -255.29 30.26
hyp7 - hyp9 -46.45 12.72
hyp7 - hyp10 -278.27 31.04
hyp8 - hyp9 208.85 31.60
hyp8 - hyp10 -22.98 9.62
hyp9 - hyp10 -231.83 28.97hyp6, 7 equivalent (no need for trt|acs) hyp8, 10 much worse (keep acs intercept random effect) hyp9 a bit worse (some evidenc for shelf effects)
go with hyp7
pred.df.species <- unique(tomato[,c("species","trt")])
hyp.species.prediction <- cbind(pred.df.species,fitted(hyp7,newdata = pred.df.species, re_formula = NA))
colnames(hyp.species.prediction)[5:6] <- c("low95","high95")
head(hyp.species.prediction)write.csv(hyp.species.prediction,file="hyp_species.csv")pl <- ggplot(hyp.species.prediction,aes(x=species,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl + ggtitle("hypocotyls per species")
ggsave("hypocotyls_species.pdf")Saving 6.77 x 4.18 in image
internode
loo(int6,int7,int8,int9,int10)Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details. LOOIC SE
int6 7371.26 65.23
int7 7419.36 64.33
int8 7635.99 66.00
int9 7433.97 62.54
int10 7644.33 64.68
int6 - int7 -48.10 17.12
int6 - int8 -264.73 38.16
int6 - int9 -62.71 19.68
int6 - int10 -273.08 38.97
int7 - int8 -216.63 27.97
int7 - int9 -14.61 9.84
int7 - int10 -224.97 29.25
int8 - int9 202.02 28.99
int8 - int10 -8.34 8.45
int9 - int10 -210.37 27.57int6 is best
int.species.prediction <- cbind(pred.df.species,fitted(int6,newdata = pred.df.species, re_formula = NA))
colnames(int.species.prediction)[5:6] <- c("low95","high95")
head(int.species.prediction)write.csv(int.species.prediction,file="int_species.csv")pl <- ggplot(int.species.prediction,aes(x=species,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl + ggtitle("internode per species")
ggsave("internodes_species.pdf")Saving 6.77 x 4.18 in image
LS0tCnRpdGxlOiAiVG9tYXRvIGZvciBDYXJsb3MiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0gCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCBjYWNoZSA9IFRSVUUsIGF1dG9kZXAgPSBUUlVFKQpgYGAKCmBgYHtyfQpsaWJyYXJ5KHJzdGFuKQpsaWJyYXJ5KHJldGhpbmtpbmcpCmxpYnJhcnkoYnJtcykKcnN0YW5fb3B0aW9ucyhhdXRvX3dyaXRlID0gVFJVRSkKb3B0aW9ucyhtYy5jb3JlcyA9IHBhcmFsbGVsOjpkZXRlY3RDb3JlcygpKQppZihmaWxlLmV4aXN0cygifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikpIGxvYWQoIn4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CnRvbWF0byA8LSByZWFkLmNzdigiVG9tYXRvUjJDU0hMLmNzdiIpCnN1bW1hcnkodG9tYXRvKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwMSA8LSBicm0oaHlwIH4gYWNzKnRydCArICh0cnR8c3BlY2llcykgKyAoMXxzaGVsZiksCiAgICAgICAgICAgIHByaW9yID0gYygKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgyMCwxMCkiLGNsYXNzPSJJbnRlcmNlcHQiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEwKSIsY2xhc3M9ImIiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoImNhdWNoeSgwLDEpIixjbGFzcz0ic2lnbWEiKQogICAgICAgICAgICApLAogICAgICAgICAgICBpdGVyPTUwMDAsCiAgICAgICAgICAgIHdhcm11cD0xNTAwLAogICAgICAgICAgICBkYXRhPXRvbWF0bykKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQpoeXAyIDwtIGJybShoeXAgfiBhY3MqdHJ0ICsgKDF8c3BlY2llcykgKyAoMXxzaGVsZiksCiAgICAgICAgICAgIHByaW9yID0gYygKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgyMCwxMCkiLGNsYXNzPSJJbnRlcmNlcHQiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEwKSIsY2xhc3M9ImIiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoImNhdWNoeSgwLDEpIixjbGFzcz0ic2lnbWEiKQogICAgICAgICAgICApLAogICAgICAgICAgICBpdGVyPTUwMDAsCiAgICAgICAgICAgIHdhcm11cD0xNTAwLAogICAgICAgICAgICBkYXRhPXRvbWF0bykKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQpoeXAzIDwtIHVwZGF0ZShoeXAxLCBmb3JtdWxhLiA9IGh5cCB+IGFjcyp0cnQgKyAoMXxzcGVjaWVzKSkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQpoeXA0IDwtIHVwZGF0ZShoeXAxLCBmb3JtdWxhLiA9IGh5cCB+IGFjcyp0cnQgKyAoMXxzaGVsZikpCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KaHlwNSA8LSB1cGRhdGUoaHlwMSwgZm9ybXVsYS4gPSBoeXAgfiBhY3MqdHJ0KQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CmxvbyhoeXAxLGh5cDIsaHlwMyxoeXA0LGh5cDUpCmBgYAoKU28gd2Ugc2hvdWxkIGdvIHdpdGggaHlwNCAobm8gc3BlY2llcyBlZmZlY3QpCmBgYHtyfQpzdW1tYXJ5KGh5cDQpCnBsb3QoaHlwNCwgYXNrPUZBTFNFKQpyYW5lZihoeXA0KQpgYGAKCmBgYHtyfQpwcmVkLmRmIDwtIHVuaXF1ZSh0b21hdG9bLGMoInNwZWNpZXMiLCJhY3MiLCJ0cnQiKV0pCmhlYWQocHJlZC5kZikKaHlwLnByZWRpY3Rpb24gPC0gY2JpbmQocHJlZC5kZixmaXR0ZWQoaHlwNCxuZXdkYXRhID0gcHJlZC5kZiwgcmVfZm9ybXVsYSA9IE5BKSkKY29sbmFtZXMoaHlwLnByZWRpY3Rpb24pWzY6N10gPC0gYygibG93OTUiLCJoaWdoOTUiKQpoZWFkKGh5cC5wcmVkaWN0aW9uKQp3cml0ZS5jc3YoaHlwLnByZWRpY3Rpb24sZmlsZT0iaHlwX2FjY2Vzc2lvbi5jc3YiKQpgYGAKCmBgYHtyfQpwbCA8LSBnZ3Bsb3QoaHlwLnByZWRpY3Rpb24sYWVzKHg9YWNzLHk9RXN0aW1hdGUseW1pbj1sb3c5NSx5bWF4PWhpZ2g5NSwgZmlsbD10cnQpKQpwbCA8LSBwbCArIGdlb21fYmFyKHBvc2l0aW9uPSJkb2RnZSIsc3RhdD0iaWRlbnRpdHkiKQpwbCA8LSBwbCArIGdlb21fZXJyb3JiYXIocG9zaXRpb249cG9zaXRpb25fZG9kZ2Uod2lkdGg9LjkpLCB3aWR0aD0wLjUpCnBsIDwtIHBsICsgZmFjZXRfd3JhcCggfiBzcGVjaWVzLCBuY29sPTIsIHNjYWxlcz0iZnJlZV94IikKcGwgKyBnZ3RpdGxlKCJoeXBvY290eWxzIHBlciBhY2Nlc3Npb24iKQpnZ3NhdmUoImh5cG9jb3R5bHNfYWNjZXNzaW9uLnBkZiIsaGVpZ2h0PTExLHdpZHRoPTguNSkKYGBgCgpgYGB7ciwgcmVzdWx0cz0naGlkZSd9CmJ5KHRvbWF0byRpbnRsZW5nLHRvbWF0byR0cnQsbWVhbikKaW50MSA8LSBicm0oaW50bGVuZyB+IGFjcyp0cnQgKyAodHJ0fHNwZWNpZXMpICsgKDF8c2hlbGYpLAogICAgICAgICAgICBwcmlvciA9IGMoCiAgICAgICAgICAgICAgc2V0X3ByaW9yKCJub3JtYWwoMTAsMTApIixjbGFzcz0iSW50ZXJjZXB0IiksCiAgICAgICAgICAgICAgc2V0X3ByaW9yKCJub3JtYWwoMCwyMCkiLGNsYXNzPSJiIiksCiAgICAgICAgICAgICAgc2V0X3ByaW9yKCJjYXVjaHkoMCwxKSIsY2xhc3M9InNpZ21hIikKICAgICAgICAgICAgKSwKICAgICAgICAgICAgaXRlcj01MDAwLAogICAgICAgICAgICB3YXJtdXA9MTUwMCwKICAgICAgICAgICAgZGF0YT10b21hdG8pCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KaW50MiA8LSB1cGRhdGUoaW50MSwgZm9ybXVsYS4gPSBpbnRsZW5nIH4gYWNzKnRydCArICgxfHNwZWNpZXMpICsgKDF8c2hlbGYpKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CmludDMgPC0gdXBkYXRlKGludDIsIGZvcm11bGEuID0gaW50bGVuZyB+IGFjcyp0cnQgKyAoMXxzcGVjaWVzKSkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQppbnQ0IDwtIHVwZGF0ZShpbnQzLCBmb3JtdWxhLiA9IGludGxlbmcgfiBhY3MqdHJ0ICsgKDF8c2hlbGYpKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CmludDUgPC0gdXBkYXRlKGludDEsIGZvcm11bGEuID0gaW50bGVuZyB+IGFjcyp0cnQpCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KbG9vKGludDEsaW50MixpbnQzLGludDQsaW50NSkKYGBgCgpTbyB3ZSBzaG91bGQgZ28gd2l0aCBpbnQ0IChubyBzcGVjaWVzIGVmZmVjdCkuICBDYXNlIGNvdWxkIGFsc28gYmUgbWFkZSBmb3IgaW50NSAobm8gcmFuZG9tIGVmZmVjdHMpCmBgYHtyfQpzdW1tYXJ5KGludDQpCnBsb3QoaW50NCwgYXNrPUZBTFNFKQpyYW5lZihpbnQ0KQpgYGAKCmBgYHtyfQppbnRsZW5nLnByZWRpY3Rpb24gPC0gY2JpbmQocHJlZC5kZixmaXR0ZWQoaW50NCxuZXdkYXRhID0gcHJlZC5kZiwgcmVfZm9ybXVsYSA9IE5BKSkKY29sbmFtZXMoaW50bGVuZy5wcmVkaWN0aW9uKVs2OjddIDwtIGMoImxvdzk1IiwiaGlnaDk1IikKaGVhZChpbnRsZW5nLnByZWRpY3Rpb24pCndyaXRlLmNzdihpbnRsZW5nLnByZWRpY3Rpb24sZmlsZT0iaW50bGVuZ19hY2Nlc3Npb24uY3N2IikKYGBgCgpgYGB7cn0KcGwgPC0gZ2dwbG90KGludGxlbmcucHJlZGljdGlvbixhZXMoeD1hY3MseT1Fc3RpbWF0ZSx5bWluPWxvdzk1LHltYXg9aGlnaDk1LCBmaWxsPXRydCkpCnBsIDwtIHBsICsgZ2VvbV9iYXIocG9zaXRpb249ImRvZGdlIixzdGF0PSJpZGVudGl0eSIpCnBsIDwtIHBsICsgZ2VvbV9lcnJvcmJhcihwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0uOSksIHdpZHRoPTAuNSkKcGwgPC0gcGwgKyBmYWNldF93cmFwKCB+IHNwZWNpZXMsIG5jb2w9Miwgc2NhbGVzPSJmcmVlX3giKQpwbCArIGdndGl0bGUoImludGVybm9kZSIpCmdnc2F2ZSgiaW50ZXJub2Rlc19hY2Nlc3Npb24ucGRmIixoZWlnaHQ9MTEsd2lkdGg9OC41KQpgYGAKCiMjIHNwZWNpZXMgbGV2ZWwgYW5hbHlzZXMKCiMjIyBoeXBvY290eWwKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwNiA8LSBicm0oaHlwIH4gc3BlY2llcyp0cnQgKyAodHJ0fGFjcykgKyAoMXxzaGVsZiksCiAgICAgICAgICAgIHByaW9yID0gYygKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgyMCwxMCkiLGNsYXNzPSJJbnRlcmNlcHQiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEwKSIsY2xhc3M9ImIiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoImNhdWNoeSgwLDEpIixjbGFzcz0ic2lnbWEiKQogICAgICAgICAgICApLAogICAgICAgICAgICBpdGVyPTUwMDAsCiAgICAgICAgICAgIHdhcm11cD0xNTAwLAogICAgICAgICAgICBkYXRhPXRvbWF0bykKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwNyA8LSB1cGRhdGUoaHlwNiwgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQgKyAoMXxhY3MpICsgKDF8c2hlbGYpICkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCgpgYGB7ciwgcmVzdWx0cz0naGlkZSd9Cmh5cDggPC0gdXBkYXRlKGh5cDcsIGZvcm11bGEuID0gLiB+IHNwZWNpZXMqdHJ0ICsgKDF8c2hlbGYpICkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwOSA8LSB1cGRhdGUoaHlwNywgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQgKyAoMXxhY3MpICkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwMTAgPC0gdXBkYXRlKGh5cDcsIGZvcm11bGEuID0gLiB+IHNwZWNpZXMqdHJ0KQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CnN1bW1hcnkoaHlwNikKYGBgCgpgYGB7cn0Kc3VtbWFyeShoeXA3KQpgYGAKCmBgYHtyfQpzdW1tYXJ5KGh5cDgpCmBgYAoKYGBge3J9CnN1bW1hcnkoaHlwOSkKYGBgCgpgYGB7cn0Kc3VtbWFyeShoeXAxMCkKYGBgCgpgYGB7cn0KbG9vKGh5cDYsaHlwNyxoeXA4LGh5cDksaHlwMTApCmBgYAoKaHlwNiwgNyBlcXVpdmFsZW50IChubyBuZWVkIGZvciB0cnR8YWNzKQpoeXA4LCAxMCBtdWNoIHdvcnNlIChrZWVwIGFjcyBpbnRlcmNlcHQgcmFuZG9tIGVmZmVjdCkKaHlwOSBhIGJpdCB3b3JzZSAoc29tZSBldmlkZW5jIGZvciBzaGVsZiBlZmZlY3RzKQoKZ28gd2l0aCBoeXA3CgpgYGB7cn0KcHJlZC5kZi5zcGVjaWVzIDwtIHVuaXF1ZSh0b21hdG9bLGMoInNwZWNpZXMiLCJ0cnQiKV0pCmh5cC5zcGVjaWVzLnByZWRpY3Rpb24gPC0gY2JpbmQocHJlZC5kZi5zcGVjaWVzLGZpdHRlZChoeXA3LG5ld2RhdGEgPSBwcmVkLmRmLnNwZWNpZXMsIHJlX2Zvcm11bGEgPSBOQSkpCmNvbG5hbWVzKGh5cC5zcGVjaWVzLnByZWRpY3Rpb24pWzU6Nl0gPC0gYygibG93OTUiLCJoaWdoOTUiKQpoZWFkKGh5cC5zcGVjaWVzLnByZWRpY3Rpb24pCndyaXRlLmNzdihoeXAuc3BlY2llcy5wcmVkaWN0aW9uLGZpbGU9Imh5cF9zcGVjaWVzLmNzdiIpCmBgYAoKYGBge3J9CnBsIDwtIGdncGxvdChoeXAuc3BlY2llcy5wcmVkaWN0aW9uLGFlcyh4PXNwZWNpZXMseT1Fc3RpbWF0ZSx5bWluPWxvdzk1LHltYXg9aGlnaDk1LCBmaWxsPXRydCkpCnBsIDwtIHBsICsgZ2VvbV9iYXIocG9zaXRpb249ImRvZGdlIixzdGF0PSJpZGVudGl0eSIpCnBsIDwtIHBsICsgZ2VvbV9lcnJvcmJhcihwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0uOSksIHdpZHRoPTAuNSkKcGwgKyBnZ3RpdGxlKCJoeXBvY290eWxzIHBlciBzcGVjaWVzIikKZ2dzYXZlKCJoeXBvY290eWxzX3NwZWNpZXMucGRmIikKYGBgCgojIyMgaW50ZXJub2RlCgpgYGB7ciwgcmVzdWx0cz0naGlkZSd9CmludDYgPC0gYnJtKGludGxlbmcgfiBzcGVjaWVzKnRydCArICh0cnR8YWNzKSArICgxfHNoZWxmKSwKICAgICAgICAgICAgcHJpb3IgPSBjKAogICAgICAgICAgICAgIHNldF9wcmlvcigibm9ybWFsKDIwLDEwKSIsY2xhc3M9IkludGVyY2VwdCIpLAogICAgICAgICAgICAgIHNldF9wcmlvcigibm9ybWFsKDAsMTApIixjbGFzcz0iYiIpLAogICAgICAgICAgICAgIHNldF9wcmlvcigiY2F1Y2h5KDAsMSkiLGNsYXNzPSJzaWdtYSIpCiAgICAgICAgICAgICksCiAgICAgICAgICAgIGl0ZXI9NTAwMCwKICAgICAgICAgICAgd2FybXVwPTE1MDAsCiAgICAgICAgICAgIGRhdGE9dG9tYXRvKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3IsIHJlc3VsdHM9J2hpZGUnfQppbnQ3IDwtIHVwZGF0ZShpbnQ2LCBmb3JtdWxhLiA9IC4gfiBzcGVjaWVzKnRydCArICgxfGFjcykgKyAoMXxzaGVsZikgKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaW50OCA8LSB1cGRhdGUoaW50NywgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQgKyAoMXxzaGVsZikgKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3IsIHJlc3VsdHM9J2hpZGUnfQppbnQ5IDwtIHVwZGF0ZShpbnQ3LCBmb3JtdWxhLiA9IC4gfiBzcGVjaWVzKnRydCArICgxfGFjcykgKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3IsIHJlc3VsdHM9J2hpZGUnfQppbnQxMCA8LSB1cGRhdGUoaW50NywgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQpCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KbG9vKGludDYsaW50NyxpbnQ4LGludDksaW50MTApCmBgYAoKaW50NiBpcyBiZXN0CgpgYGB7cn0KaW50LnNwZWNpZXMucHJlZGljdGlvbiA8LSBjYmluZChwcmVkLmRmLnNwZWNpZXMsZml0dGVkKGludDYsbmV3ZGF0YSA9IHByZWQuZGYuc3BlY2llcywgcmVfZm9ybXVsYSA9IE5BKSkKY29sbmFtZXMoaW50LnNwZWNpZXMucHJlZGljdGlvbilbNTo2XSA8LSBjKCJsb3c5NSIsImhpZ2g5NSIpCmhlYWQoaW50LnNwZWNpZXMucHJlZGljdGlvbikKd3JpdGUuY3N2KGludC5zcGVjaWVzLnByZWRpY3Rpb24sZmlsZT0iaW50X3NwZWNpZXMuY3N2IikKYGBgCgpgYGB7cn0KcGwgPC0gZ2dwbG90KGludC5zcGVjaWVzLnByZWRpY3Rpb24sYWVzKHg9c3BlY2llcyx5PUVzdGltYXRlLHltaW49bG93OTUseW1heD1oaWdoOTUsIGZpbGw9dHJ0KSkKcGwgPC0gcGwgKyBnZW9tX2Jhcihwb3NpdGlvbj0iZG9kZ2UiLHN0YXQ9ImlkZW50aXR5IikKcGwgPC0gcGwgKyBnZW9tX2Vycm9yYmFyKHBvc2l0aW9uPXBvc2l0aW9uX2RvZGdlKHdpZHRoPS45KSwgd2lkdGg9MC41KQpwbCArIGdndGl0bGUoImludGVybm9kZSBwZXIgc3BlY2llcyIpCmdnc2F2ZSgiaW50ZXJub2Rlc19zcGVjaWVzLnBkZiIpCmBgYA==